home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / demo3.i < prev    next >
Text File  |  1995-07-26  |  4KB  |  141 lines
  1. /*
  2.    DEMO3.I
  3.    Chaotic pendulum demo
  4.  
  5.    $Id$
  6.  */
  7. /*    Copyright (c) 1995.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. func demo3(time_limit, kick=, damp_time=)
  11. /* DOCUMENT demo3
  12.      Solve Lagrange's equations for a famous chaotic pendulum (as on
  13.      exhibit at the San Francisco Exploratorium).  Run a movie of the
  14.      result.  By default, the movie runs for 60 seconds, but if you
  15.      supply an argument to the demo3 function, it will run for that
  16.      many seconds.
  17.  
  18.      The kick= keyword may be used to adjsut the initial amplitude;
  19.      kick= 1.2 is the default.
  20.  
  21.      You may also wish to supply a damp_time= keyword to see the effect
  22.      of an ad hoc damping term.  damp_time=100 is nice.
  23.  */
  24. {
  25.   require, "rkutta.i";
  26.   a= 2.0;
  27.   b= 3.0;
  28.   c= 1.5;
  29.   gravity= 1.0;
  30.   widget_setup, a, b, c, 1.00, 0.67, gravity;
  31.  
  32.   dt= 0.05*sqrt(min(a,b,c)/gravity);
  33.  
  34.   if (is_void(damp_time)) damp_time= 1.e30;
  35.   if (is_void(kick)) kick= 1.2;
  36.  
  37.   /* roll dice to get initial velocities */
  38.   thetadot= kick;
  39.   gammadot= 0.1*random();
  40.   qdotq= [[thetadot,0.,0.,gammadot], [0.,0.,0.,0.]];
  41.  
  42.   window, wait=1;
  43.   s= 1.1*(b+c);
  44.   limits, -s,s,-s,s;
  45.   display, qdotq;  /* without this, won't get axis numbers */
  46.  
  47.   require, "movie.i";
  48.   movie, draw_frame, time_limit, 0.02;
  49. }
  50.  
  51. func draw_frame(i)
  52. {
  53.   display, qdotq;
  54.   qdotq= rkdumb(lagrange, qdotq,0., dt, 1, nosave=1);
  55.   return 1;  /* just go until time limit expires */
  56. }
  57.  
  58. func widget_setup(a, b, c, m_tee, m_bar, gravity)
  59. {
  60.   /* set up the values used by lagrange */
  61.   extern Tdiag, Ta, Tb, Tc, Ba, Bb, Bg, Dg;
  62.  
  63.   B= 0.5*m_bar*c^2;
  64.   C= 0.5*B*c;
  65.   A= m_tee*(2.*a^3+b^3)/3. + m_bar*c*(2.*a^2+b^2) + C;
  66.   Ba= B*a;
  67.   Bb= B*b;
  68.  
  69.   Dg= (0.5*m_tee*b^2 + m_bar*b*c)*gravity;
  70.   Bg= B*gravity;
  71.  
  72.   Tdiag= [[A,0.,0.,0.], [0.,C,0.,0.], [0.,0.,C,0.], [0.,0.,0.,C]];
  73.   Ta= [[0.,Ba,0.,0.], [Ba,0.,0.,0.], [0.,0.,0.,0.], [0.,0.,0.,0.]];
  74.   Tb= [[0.,0.,Ba,0.], [0.,0.,0.,0.], [Ba,0.,0.,0.], [0.,0.,0.,0.]];
  75.   Tc= [[0.,0.,0.,Bb], [0.,0.,0.,0.], [0.,0.,0.,0.], [Bb,0.,0.,0.]];
  76. }
  77.  
  78. func lagrange(qdotq, time/*unused*/)
  79. {
  80.   /* compute the time derivatives of [qdot,q]
  81.      Lagrange's equations d(dL/dqdot)/dt - dL/dq = 0
  82.      turn out to be T*qdotdot+L=0, with T and L computed here.  */
  83.   theta= qdotq(1,2);
  84.   alpha= qdotq(2,2);
  85.   beta= qdotq(3,2);
  86.   gamma= qdotq(4,2);
  87.   qdot= qdotq(,1);
  88.   t2= qdot(1)^2;
  89.   a2= qdot(2)^2;
  90.   b2= qdot(3)^2;
  91.   g2= qdot(4)^2;
  92.  
  93.   amt= alpha-theta;
  94.   bpt= beta+theta;
  95.   gmt= gamma-theta;
  96.  
  97.   T= Tdiag + Ta*sin(amt) + Tb*sin(bpt) - Tc*cos(gmt);
  98.   camt= cos(amt);
  99.   cbpt= cos(bpt);
  100.   sgmt= sin(gmt);
  101.   L= [Dg*sin(theta) + Ba*(a2*camt+b2*cbpt)+Bb*g2*sgmt,
  102.       Bg*sin(alpha) - Ba*t2*camt,
  103.       Bg*sin(beta) + Ba*t2*cbpt,
  104.       Bg*sin(gamma) - Bb*t2*sgmt];
  105.   qdotdot= LUsolve(T, -L);
  106.  
  107.   qdotdot(,1)-= qdot/damp_time;
  108.  
  109.   return [qdotdot, qdot];
  110. }
  111.  
  112. func display(qdotq)
  113. {
  114.   /* draw the widget given its current state */
  115.   theta= qdotq(1,2);
  116.   alpha= qdotq(2,2);
  117.   beta= qdotq(3,2);
  118.   gamma= qdotq(4,2);
  119.  
  120.   st= sin(theta);
  121.   ct= cos(theta);
  122.   p1x= -a*ct;
  123.   p1y= -a*st;
  124.   p2x= -p1x;
  125.   p2y= -p1y;
  126.   p3x= b*st;
  127.   p3y= -b*ct;
  128.   q1x= p1x + c*sin(alpha);
  129.   q1y= p1y - c*cos(alpha);
  130.   q2x= p2x + c*sin(beta);
  131.   q2y= p2y - c*cos(beta);
  132.   q3x= p3x + c*sin(gamma);
  133.   q3y= p3y - c*cos(gamma);
  134.  
  135.   plg, [p1y,p2y], [p1x,p2x], marks=0,type=1,width=36,color="red";
  136.   plg, [0.,p3y], [0.,p3x], marks=0,type=1,width=36,color="red";
  137.   plg, [p1y,q1y], [p1x,q1x], marks=0,type=1,width=24,color="green";
  138.   plg, [p2y,q2y], [p2x,q2x], marks=0,type=1,width=24,color="green";
  139.   plg, [p3y,q3y], [p3x,q3x], marks=0,type=1,width=24,color="green";
  140. }
  141.